function out = AR1_regimeSwitching(yields,maturities,econAct,thresholdLevel)

regime_rec     = real(econAct) < thresholdLevel;  %a recression
[T,nMat]       = size(yields);
bettaHat       = zeros(nMat,4);
tstatHat       = zeros(nMat,4);
seBettaHat     = zeros(nMat,4);
R2             = zeros(nMat,1);
for i=1:nMat
    Y = yields(2:end,i);
    X = [ones(T-1,1) regime_rec(1:end-1,1).*ones(T-1,1) yields(1:end-1,i)  regime_rec(1:end-1,1).*yields(1:end-1,i)];
    res = nwest(Y,X,3+3+6);
    bettaHat(i,:)   = res.beta(1:4)';
    tstatHat(i,:)   = res.tstat(1:4)';
    seBettaHat(i,:) = res.beta(1:4)./res.tstat(1:4);
    R2(i,1)         = res.rbar;
end

% Output
out.bettaHat   = bettaHat;
out.seBettaHat = seBettaHat;
out.tstatHat   = tstatHat;
out.maturities = maturities;
out.R2         = R2;

end

